testP = rand(100, 2);

for i = 1 : 100
    F1 = evaluateInteriorDeformationGradient(testP(i,1), testP(i,2), domain2D, theta);
    F2 = getDeformationGradientFD(testP(i,1), testP(i,2), domain2D, theta);
    
    F1 + F1' - F2
end